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Abstract 

The dual-phase lagging (DPL) model has been considered as one of the most promising theoretical approaches to 
generalize the classical Fourier law for heat conduction involving short time and space scales. Its applicability, 
potential, equivalences, and possible drawbacks have been discussed in the current literature. In this study, the 
implications of solving the exact DPL model of heat conduction in a three-dimensional bounded domain solution 
are explored. Based on the principle of causality, it is shown that the temperature gradient must be always the 
cause and the heat flux must be the effect in the process of heat transfer under the dual-phase model. This fact 
establishes explicitly that the single- and DPL models with different physical origins are mathematically equivalent. 
In addition, taking into account the properties of the Lambert W function and by requiring that the temperature 
remains stable, in such a way that it does not go to infinity when the time increases, it is shown that the DPL 
model in its exact form cannot provide a general description of the heat conduction phenomena. 



Introduction 

Nanoscale heat transfer involves a highly complex pro- 
cess, as has been witnessed in the last years in which 
remarkable novel phenomena related to very short time 
and spatial scales, such as enhancement of thermal con- 
ductivity in nanofluids, granular materials, thin layers, 
and composite systems among others, have been 
reported [1-5]. The traditional approach to deal with 
these phenomena has been to use the Fourier heat trans- 
fer equation. This methodology has proven to be exten- 
sively useful in the analysis of heat transport in a great 
variety of physical systems, however, when applied to 
highly heterogeneous systems or when the time and 
space scale are very short, they show serious inconsisten- 
cies [6,7]. In order to understand the nanoscale heat 
transfer, a great diversity of novel theoretical approaches 
have been developed [3,5,7,8]. In particular, when analyz- 
ing two-phase systems, one of the simplest heat conduc- 
tion models that considers the microstructure is known 
as the two-equation model [9,10], which has been devel- 
oped writing the Fourier law of heat conduction [11] for 
each phase and performing a volume averaging proce- 
dure [9]. This model takes into account the porosity of 
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the component phases as well as their interface effects by 
means of two coefficients [12]. Besides, it has been 
shown that the two-equation model is equivalent to the 
one-equation model known as the dual-phase lagging 
(DPL) model, in which the microstructural effects 
are taken into account by means of two time delays 
[3,10,13-15]. DPL model have been proposed to sur- 
mount the well-known drawbacks of the Fourier law and 
the Cattaneo equation of heat conduction [7], and estab- 
lishes that either the temperature gradient may precede 
the heat flux or the heat flux may precede the tempera- 
ture gradient. Mathematically, this is written in the form 



Cj(x, t + Tq) = — feVT(x, t + Tj), 



(1) 



where x is the position vector, t is the time, 
cj [W- m~ 2 ] is the heat flux vector, T[K] is the absolute 
temperature, /r[W.nT 1 .K~ 1 ] is the thermal conductivity, 
x q is the phase lag of the heat flux, and x T is the phase 
lag of the temperature gradient. For the case of x q >Xr> 
the heat flux (effect) established across the material is a 
result of the temperature gradient (cause); while for 
i q <%T> the heat flux (cause) induces the temperature gra- 
dient (effect). Notice that when x q = x T , the response 
between the temperature gradient and the heat flux is 
instantaneous and Equation 1 reduces to Fourier law 
except for a trivial shift in the time scale. In addition, 
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note that for x T = 0; the DPL model reduces to the sin- 
gle-phase lagging (SPL) model [3], The time delay x q is 
interpreted as the relaxation time due to the fast-transi- 
ent effects of thermal inertia, while the phase lag x T 
represents the time required for the thermal activation 
in micro-scale [3]. For the case of composite materials, 
the phase lag x q takes into account the time delay due 
to contact thermal resistance among the particles, while 
x T is interpreted as the time required to establish the 
temperature gradient through the particles [12,16]. The 
lagging behavior in the transient process is caused by 
the finite time required for the microscopic interactions 
to take place. This time of response has been claimed to 
be in the range of a few nanoseconds in metals and up 
to the order of several seconds in granular matter [3]. In 
this last case, due to the low-conducting pores among 
the grains and their interface thermal resistance. 

The thermal conductivity is an intrinsic property of 
each material which measures its ability for the transfer 
of heat and is determined by the kinetic properties of 
the energy carriers and the material microstructure 
[6,17]. Under the framework of Boltzmann kinetic the- 
ory [3,6], it can be shown that the thermal conductivity 
is directly proportional to the group velocity and mean 
free path of the energy carriers (electrons and phonons). 
These parameters depend strongly on the material tem- 
perature, due to the multiple scattering processes 
involved among energy carriers and defects, such as 
impurities, dislocations, and grain boundaries, [6,18]. 
Thus, in general; thermal conductivity exhibits compli- 
cated temperature dependence. However, in many cases 
of practical interest, the thermal conductivity can be 
considered independent of the temperature for a consid- 
erable range of operating temperatures [3,6,11]. Based 
on this fact and to keep our mathematical approach 
tractable, we assume that thermal conductivity is a tem- 
perature-independent parameter. 

Phase lags represent the time parameters required by 
the material to start up the heat flux and temperature 
gradient, after a thermal excitation has been imposed; 
larger phase lags are expected in material with smaller 
thermal conductivities, as is the case of granular matter 
[3]. Materials, in which the temperature gradient phase 
lag dominates, show a strong attenuation of the neat 
heat flux. In this case, the behavior is dominated by 
parabolic terms of the heat transport equation. In con- 
trast, materials in which the heat flux phase lag is domi- 
nant show a slight attenuation of the heat flux, implying 
that a hyperbolic Cattaneo-Vernotte heat propagation is 
present. For a further discussion of the relationship 
between thermal conductivity and phase lags, Tzou's 
book [3] is recommended. 

It is convenient to take into account that the heat flux 
and temperature gradient shown in Equation 1 are the 



local responses within the medium. They must not be 
confused with the global quantities specified in the 
boundary conditions. When a heat flux (as a laser 
source) is applied to the boundary of a solid medium, 
the temperature gradient established within the medium 
can still precede the heat flux. The application of the 
heat flux at the boundary does not guarantee the prece- 
dence of the heat flux vector to the temperature gradi- 
ent at all. In fact, whether the heat flux vector precedes 
the temperature gradient or not depends on the com- 
bined effects of the thermal loading and thermal proper- 
ties of the materials, as was explained by Tzou [3]. In 
this way, the DPL model should provide a comprehen- 
sive treatment of the heterogeneous nature of composite 
media [3,13]. 

It has been shown that under the DPL model and in 
absence of internal heat sources, the temperature satis- 
fies the following differential-difference equation 
[19-22]: 

V 2 T(S, t -r)-i^M = 0, (2) 
a dt 

where atm^s' 1 ] is the thermal diffusivity of the med- 
ium, and x = x q -x T is the difference of the phase lags. 
Equation 2 shows explicitly that the DPL and SPL mod- 
els, both in their exact form, are entirely equivalent, 
when x> 0(x q -x T )[19]. 

The solutions of Equation 2 for some geometries have 
been explored [19-22]. In the time domain, Jordan et al. 
[19] and Quintanilla and Jordan [22] have shown that 
the SPL model, in its exact form, can lead to instabilities 
with respect to specific initial values. Additionally, in the 
frequency domain, using a modulated heat source, 
Ordonez-Miranda and Alvarado-Gil [21] have shown 
that the if the DPL model is valid, its applicability must 
be restricted to frequency-interval strips, which are 
determined only by the difference of the time delays x = 
x q -x T . These studies have pointed out that the usefulness 
of the Cattaneo-Vernotte and DPL exact models is 
limited. 

In this study, by means of the method of separation of 
variables, the solution of Equation 2 is obtained in a 
bounded domain. It is shown that, for any kind of 
homogeneous boundary conditions, its solutions go to 
infinity in the long time domain. This explosive charac- 
teristic of the temperature predicted by Equation 2 indi- 
cates that the DPL model, in its exact form, cannot be 
considered as a valid model of heat conduction. 

Mathematical formulation and solutions 

The general solution of Equation 2 in a three-dimen- 
sional closed region of finite volume V and boundary 
surface dV is going to be obtained in this section. The 
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initial-boundary value problem to be solved can be writ- 
ten as follows: 

V 2 T(x, £ - t) - - 9T ^' - = 0, (x, t) e V x (0, +oo); (3a) 
a dt 



aT(x, t) + bVT{x, t) ■ n = 0, (x, t) e dV x (0, +oo); (3b) 



T(x, t) = T 0 {x, t), {x, t) € V x [-t, 0]; 



(3c) 



where a and & are two constants and n is a unit nor- 
mal vector pointing outward of the boundary surface 
dV. Note that the boundary conditions in Equation 3a 
imply the specification of the temperature and heat flux 
at dV and they reduce to the Dirichlet (Neumann) pro- 
blem for b = 0 {a = 0) [5]. On the other hand, the initial 
condition is specified in the pre-interval [-x,0] to define 
the time derivative of the temperature in the interval [0, 
x]. This is a common characteristic of the delay differen- 
tial equations, as Equation 3a [23]. In many common 
situations the initial history function Tn(x,t) may be 
considered as a constant. 

According to the method of separation of variables, a 
solution of the form 



T(x,t) = f(x)p(t), 



(4) 



is proposed. After inserting Equation 4 into Equations 
3a, b, it is obtained that 



V ir m {x) +k m ir n (x) = 0, 



ai/r m (x) + bVi/ m (x) ■ h = 0, 



dpmjt) 

dt 



+ Uk m p m {t- T) = 0, 



(5a) 



(5b) 



(5c) 



where the integer subscript m = 1,2,3,... has been 
inserted in view that Equations 5a, b defined an 
eigenvalue (Sturm-Liouville) problem [5], and X m is 
the eigenvalue associated with the eigenfunction y/ m . 
As an example, in the case of one-dimensional heat 
conduction across a finite region 0 <x<l, nine possi- 
ble combinations of the boundary conditions given 
by Equation 5b can be found [5]. One of these com- 
binations occurs when both surfaces x = 0 and x = I 
are insu lated (d^/dbc| x=0 = di/r/dx| j=; = 0). After 
applying these particular boundary conditions to the 
solution of Equation 5a, it is found that its eigenva- 
lues are determined by \ m = (mn/l) 2 - Similar results 
can be obtained for the other combinations of 
boundary conditions as well as for more complex 
geometries [5]. In general, all the eigenvalues are 
real and positive, and they go to infinity when 



m—>°°[5]. In this way, by the principle of superposi- 
tion, the general solution of Equation 3a-c can be 
written as 



T(x,i) = ^2if m {x)p m {t), 



(6) 



where Equation 5c can be solved assuming that P m {t) 
= exp{st) is its solution for some value of 5. This pro- 
vides the relationship 



s + aX m e 



0, 



(7) 



whose solutions can be expressed in a closed form by 
means of the Lambert W function as follows [24]: 



s m , r r = W r {-arX m ), 



(8) 



where r = 0,± 1,± 2,... indicates a specific branch of 
the complex-valued function W r . For y^-e' 1 , all the 
branches of W r (y) are different; while for y = -e , the 
branches W.\{y) = W 0 (y) = -1 and the others have dif- 
ferent values among them. In this way, the general solu- 
tion of Equation 5c is given by 

+00 

Pm(t)= C m , r exp[W r (-ar* m )t/T], m j M, (9a) 



,,(() = (D„,, 0 +D,„,_if)exp(-(/r) + J] D,„, r exp [W, .{-e ')I/r], m = M fcfo^ 



where axk M = e 1 and the constants C m>r and D m r can 
be determined by expanding Equation 3c in terms of 
the orthogonal set of eigenfunctions {y/ m } as follows: 



T 0 (x, t) = ^b m (i)f m (x). 

m=\ 

In this way, for -x<t< 0 
Pm{t) = b m (i), 



(10) 



(11) 



is satisfied. However, in practice the determination of 
the coefficients C m , r and D m r by means of Equation 11 
may be complicated. This can be avoided by solving 
Equation 5c using the Laplace transform method. After 
taking the Laplace transform of Equation 5c, and using 
Equation 11, it is obtained that in the Laplace domain, 
the function P m (s)=L[P m (t)] is given by 



b m (0) - aX m B m {s)e 
s + ak m e~ sz 



(12) 



where B m (s)=L[b m (t)] for the time domain -x<t< 0. 
Using the complex inversion formula of the Laplace 
transform [5], it is obtained that 
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1 {t) = J2 R [ P m{s)e st ,S = Sr, m ], 



(13) 



where R[] stands for the residue of its argument. 
Given that the poles of Equation 12 are determined by 
equating to zero its denominator, these poles s r>m are 
determined by Equation 8. Note that all the poles are 
simple if axX m *e x , and there is a double pole for axX m 
= e' 1 , at r = -1,0. In this way, after calculating the resi- 
dues involved in Equation 13 and comparing Equations 
9a, b with Equation 13 it is found that 



b m (0) + S m ,,-B m (s m , r ) 



1 +S„ 



(14a) 



fe m (0) + r- 1 W,-(-e- 1 )B m [r- 1 W^-e' 1 )] 
1 + Wri-e- 1 ) 



(14b) 



D m , 0 = - [MO) + 2r- 1 B,„ (-r- 1 ) - 3r- 2 B' m (-I" 1 )] ,(14c) 



D m ,_! = It' 1 [b m {0) - x- l B m (-r^ 1 )] , 



(14d) 



where the parameters s r>m are given by Equation 8 
and the prime (') on B m indicates derivative with 
respect to its argument. For the particular case in 
which the initial history function does not depend on 
time, the coefficient b m = constant =b 0 and Equations 
14a-d reduce to 



Cm 



D„ 



-b 0 aX„ 



b 0 e 



(15a) 



WrC-e-^fl + WrC-e- 1 )]' 



8e- 



2<T 1 r- 1 f> 0 , 



(15c) 



(15d) 



which agree with the previous results of Jordan et al. 
[19]. It is interesting to note that by requiring that P m 
(0) = b 0 in Equation 9a, the following property of the 
Lambert W function is obtained 



E 



i 



W r (y) [1 + W r (y)] 



(16) 



where y = -axX m . Using appropriate software, Equa- 
tion 16 can be verified to be valid not only for the roots 
of Equation 7, but also for any value of y. 



Analysis of the results 

In this section, the time-dependent part of the tempera- 
ture is going to be analyzed in two key points, as follows: 

• According to Equation 5c, the temporal rate of 
change of P m (t) (and therefore of the temperature) is 
determined by its value at the past (future), if x> 0 (t< 0). 
Based on the principle of causality, the future cannot 
determine the past, and therefore the DPL model in its 
exact form (Equation 1) must take into account the con- 
straint t = x q -x T > 0. In this way, the DPL and SPL models 
are fully equivalent between them [3,5]. This fact is in 
strong contrast to the values of the phase lags, reported 
by Tzou [3]. By expanding both sides of Equation 1 in a 
Taylor series and considering a first-order approximation 
in the phase lags, this author found that x T = 100 x q for 
metals. This discrepancy with the causality principle indi- 
cates that the predictions of the DPL model in its approx- 
imate and exact forms may be remarkably different. This 
fact reveals that the small-phase lags can have great 
effects, as it has been shown in the theory of delayed dif- 
ferential equations [23]. 

♦ Based on Equation 9a and taking into account that 
the principle of causality demands that x> 0, as has been 
discussed in above, it can be observed that the tempera- 
ture remains stable (finite) for large values of time, if 
the following condition is satisfied 



Re [W, (-arA m )] < 0, Vr e Z; 



(17) 



where Re[] stands for the real part of its argument. 
For y = Til 2, Figure 1 shows that the larger real parts of 
W r (y) are given when r = -1,0. In general, after a graphi- 
cal analysis of the Lambert W function, it can be con- 
cluded that max{Re[W r (y)]} = W 0 {y)Xy e 0t[24]. 
Based on this result, Equation 17 can be replaced by 



(15b) Re[W Q (-axX m )] < 0. 



(18) 



Given that Re [W 0 (y < > 0„ 

Re[W 0 {-jr/2 < y < 0)] < 0 and 
Re[W 0 {y = -jr/2)] = 0 (see Fi gure 1), the inequality 
(18) is satisfied if and only if 



7T 

aTk m < —, 
2 



Vm = 1,2,3,... 



(19) 



which represent the stability condition of the tempera- 
ture for long times. 

Taking into account that X m — »°° for m— >°°, it can be 
observed that the condition (19) cannot be satisfied for 
arbitrarily large values of m. The only way to solve this 
would be by imposing that m<m max , in such a way that 
aTk m = 7T 1 2, however, under this restriction on the 

'"max / ' ' 

values of m, the initial condition could not be satisfied 
(Equation 10). In this way, it is concluded that the DPL 
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Figure 1 Distribution of the imaginary values of W r (y) with respect to its real values, at y = n/2. 



model in its exact form establishes that the temperature 
increases without limit when the time grows, which is phy- 
sically unacceptable. This divergent behavior of the tem- 
perature, in the DPL model at long times, is the direct 
consequence of having introduced the phase lags. Even 
though the effects of these parameters are obviously very 
important for short time scales, according to our results 
(see Equations 1 and 9a, b), the assumption of talcing them 
as different from zero implies non-physical behavior at 
large time scales. Therefore, the DPL model, in its exact 
form, cannot be a valid formalism for heat conduction 
analysis in the complete time scale. It is expected that the 
correct model of heat conduction at both short and large 
scales could be derived from the Boltzmann transport 
equation under the relaxation time approximation [6] . 

Conclusions 

By combining the methods of separation of variables 
and the Laplace transform, the exact solution of the 
DPL model of heat conduction in a three-dimensional 



bounded domain has been obtained and analyzed. 
According to the principle of causality, it has been 
shown that the temperature gradient must precede the 
heat flux. In addition, based on the properties of the 
Lambert W function, it has been shown that the DPL 
model predicts that the temperature increases without 
limit when the time goes to infinity. This unrealistic 
prediction indicates that the DPL model, in its exact 
form, does not provide a general description of the heat 
conduction phenomena for all time scales as had been 
previously proposed. 
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